#!/usr/bin/perl -w
use strict;
use Time::HiRes q/gettimeofday/;
use POSIX q/strftime/;


my $samtools="/bluearc/data/schatz/software/bin/samtools";
#my $FA_path="/bluearc/data/schatz/hadoop/chr22_org/chr22.fa";
my $FA_path="chr22.fa";
my $bam_name = "out.bam";
my $sam_name = "out.sam";
my $sorted_bam_name ="out.prefix";
my $pileup_fname = "pileup";
my $bam_h = "bam/";

my ($seconds, $microseconds) = gettimeofday;
my $counter = strftime( q/%Y%m%dT%H%M%S/, localtime($seconds)),
	     sprintf("%05d", $microseconds/10);
my $outsampath = $counter.$sam_name;

my %hash;
my $file_counter = 0;

while(<>){
        chomp;
	my $line = $_;
	
	my $fh;
	$line =~ /(\@\S+)(\s+?)(\S+?)(\s+?)(\S+?)(\s+?)(.*)/;	

	#print "###\n\$1: $1\n\$2: $2\n\$3: $3\n\$4: $4\n\$5: $5\n\$6: $6\n\$7: $7\n###\n";	
	my $key = $1.$2.$3.$4.$5;
	my $cur_outsampath;
	if(!exists $hash{$key}){
		$file_counter++;
		$hash{$key} = $file_counter;
		$cur_outsampath = $hash{$key}.$outsampath;
		open $fh,">>".$cur_outsampath or die "could not open at $!";

       		print $fh $key;
	        close $fh; 
	}
	$cur_outsampath = $hash{$key}.$outsampath;
	open $fh,">>".$cur_outsampath or die "could not open at $!";

       	print $fh "\n$7";
        close $fh; 
}

foreach my $key(keys %hash){
	my $cur_outsampath = $hash{$key}.$outsampath;
	my $fh;
        open $fh,">>".$cur_outsampath or die "could not open at $!";

       	print $fh "\n";
        close $fh; 
	my $outbampath = $cur_outsampath.$counter.$bam_name;
        system("$samtools view -bSu $cur_outsampath > $outbampath") == 0 or die "could not execute at $!";

	#print STDERR "\$outbampath: $outbampath\n";
        my $sortedbampath = $cur_outsampath.$counter.$sorted_bam_name;

        system("$samtools sort $outbampath $sortedbampath") == 0 or die "could not execute at $!";

        $sortedbampath .= ".bam";
        system("$samtools index $sortedbampath") == 0 or die "could not execute at $!";

	$pileup_fname = $cur_outsampath.$counter.$pileup_fname;
	system("$samtools mpileup -f $FA_path $sortedbampath > $pileup_fname") == 0 or die "could not execute at $!";
	
	open pileup_fh,$pileup_fname;
	while(<pileup_fh>){
		print $_;
	}
	close pileup_fh;
}
	



print STDERR "END_OF wrap_pileup.pl \n";
